home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / vvector.cc < prev    next >
Encoding:
Text File  |  1995-12-20  |  15.3 KB  |  529 lines  |  [TEXT/ALFA]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *            Verify Operations on Vectors
  6.  *
  7.  * $Id: vvector.cc,v 3.1 1995/02/01 16:37:40 oleg Exp oleg $
  8.  *
  9.  ************************************************************************
  10.  */
  11.  
  12. #include "LinAlg.h"
  13. #include <math.h>
  14. #include "builtin.h"
  15. #include <iostream.h>
  16.  
  17. /*
  18.  *------------------------------------------------------------------------
  19.  *            Service validation functions
  20.  */
  21. static void verify_vector_identity(const Vector& v1, const Vector& v2)
  22. { verify_matrix_identity(v1,v2); }
  23.  
  24. /*
  25.  *------------------------------------------------------------------------
  26.  *       Test allocation functions and compatibility check
  27.  */
  28.  
  29. static void test_allocation(void)
  30. {
  31.  
  32.   cout << "\n\n---> Test allocation and compatibility check\n";
  33.  
  34.   Vector v1(20);
  35.   Vector v2(1,20);
  36.   Vector v3(0,19);
  37.   Vector v4(v1);
  38.   cout << "The following vector have been allocated\n";
  39.   v1.info(), v2.info(), v3.info(), v4.info();
  40.  
  41.   cout << "\nStatus information reported for vector v3:\n";
  42.   cout << "  Lower bound ... " << v3.q_lwb() << "\n";
  43.   cout << "  Upper bound ... " << v3.q_upb() << "\n";
  44.   cout << "  No. of elements " << v3.q_no_elems() << "\n";
  45.   cout << "  Name " << v3.q_name() << "\n";
  46.  
  47.   cout << "\nCheck vectors 1 & 2 for compatibility\n";
  48.   are_compatible(v1,v2);
  49.  
  50.   cout << "Check vectors 1 & 4 for compatibility\n";
  51.   are_compatible(v1,v4);
  52.  
  53. #if 0
  54.   cout << "Check vectors 1 & 3 for compatibility\n";
  55.   are_compatible(v1,v3);
  56. #endif
  57.  
  58.   cout << "v2 has to be compatible with v3 after resizing to v3\n";
  59.   v2.resize_to(v3);
  60.   are_compatible(v2,v3);
  61.  
  62.   Vector v5(v1.q_upb()+5); v5.set_name("Vector v5");
  63.   cout << "v1 has to be compatible with v5 after resizing to v5.upb\n";
  64.   v5.info(); cout << endl;
  65.   v1.resize_to(v5.q_no_elems());
  66.   are_compatible(v1,v5);
  67.  
  68.   {
  69.     cout << "Check that shrinking does not change remaining elements" << endl;
  70.     Vector vb(-1,20);
  71.     register int i;
  72.     for(i=vb.q_lwb(); i<=vb.q_upb(); i++)
  73.       vb(i) = i+M_PI;
  74.     Vector v = vb;
  75.     assert( v == vb );
  76.     assert( v != 0 );
  77.     v.resize_to(0,10);
  78.     for(i=v.q_lwb(); i<=v.q_upb(); i++)
  79.       assert( v(i) == vb(i-v.q_lwb()+vb.q_lwb()) );
  80.     cout << "Check that expansion expands by zeros" << endl;
  81.     const int old_nelems = v.q_upb() - v.q_lwb() + 1;
  82.     v.resize_to(vb);
  83.     assert( !(v == vb) );
  84.     for(i=v.q_lwb(); i<v.q_lwb()+old_nelems; i++)
  85.       assert( v(i) == vb(i) );
  86.     for(; i<=v.q_upb(); i++)
  87.       assert( v(i) == 0 );
  88.   }
  89.   cout << "\nDone\n";
  90. }
  91.  
  92. /*
  93.  *------------------------------------------------------------------------
  94.  *             Test uniform element operations
  95.  */
  96.  
  97. static void test_element_op(const int vsize)
  98. {
  99.   const double pattern = M_PI;
  100.   register int i;
  101.  
  102.   cout << "\n\n---> Test operations that treat each element uniformly\n";
  103.  
  104.   Vector v(-1,vsize-2);
  105.   Vector v1(v);
  106.  
  107.   cout << "Writing zeros to v...\n";
  108.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  109.     v(i) = 0;
  110.   verify_element_value(v,0);
  111.  
  112.   cout << "Clearing v1 ...\n";
  113.   v1.clear();
  114.   verify_element_value(v1,0);
  115.  
  116.   cout << "Comparing v1 with 0 ...\n";
  117.   assure(v1 == 0, "v1 is not zero!");
  118.  
  119.   cout << "Writing a pattern " << pattern << " by assigning to v(i)...\n";
  120.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  121.     v(i) = pattern;
  122.   verify_element_value(v,pattern);
  123.  
  124.   cout << "Writing the pattern by assigning to v1 as a whole ...\n";
  125.   v1 = pattern;
  126.   verify_element_value(v1,pattern);
  127.  
  128.   cout << "Comparing v and v1 ...\n";
  129.   assure(v == v1, "v and v1 are unexpectedly different!");
  130.   cout << "Comparing (v=0) and v1 ...\n";
  131.   assert(!(v.clear() == v1));
  132.  
  133.   cout << "\nClear v and add the pattern\n";
  134.   v.clear();
  135.   v += pattern;
  136.   verify_element_value(v,pattern);
  137.   cout << "   add the doubled pattern with the negative sign\n";
  138.   v += -2*pattern;
  139.   verify_element_value(v,-pattern);
  140.   cout << "   subtract the trippled pattern with the negative sign\n";
  141.   v -= -3*pattern;
  142.   verify_element_value(v,2*pattern);
  143.  
  144.   cout << "\nVerify comparison operations\n";
  145.   v = pattern;
  146.   assert( v == pattern && !(v != pattern) && v >= pattern && v <= pattern );
  147.   assert( v > 0 && v >= 0 );
  148.   assert( v > -pattern && v >= -pattern );
  149.   assert( v < pattern+1 && v <= pattern+1 );
  150.   v(v.q_upb()) += 1;
  151.   assert( !(v==pattern) && !(v != pattern) && v != pattern-1 );
  152.   assert( v >= pattern && !(v > pattern) && !(v >= pattern+1) );
  153.   assert( v <= (REAL)pattern+1 && !(v < pattern+1) && !(v <= pattern) );
  154.  
  155.   cout << "\nAssign 2*pattern to v by repeating additions\n";
  156.   v = 0; v += pattern; v += pattern;
  157.   cout << "Assign 2*pattern to v1 by multiplying by two \n";
  158.   v1 = pattern; v1 *= 2;
  159.   verify_element_value(v1,2*pattern);
  160.   assert( v == v1 );
  161.   cout << "Multiply v1 by one half returning it to the 1*pattern\n";
  162.   v1 *= 1/2.;
  163.   verify_element_value(v1,pattern);
  164.  
  165.   cout << "\nAssign -pattern to v and v1\n";
  166.   v.clear(); v -= pattern; v1 = -pattern;
  167.   verify_element_value(v,-pattern);
  168.   assert( v == v1 );
  169.   cout << "v = sqrt(sqr(v)); v1 = abs(v1); Now v and v1 have to be the same\n";
  170.   v.sqr();
  171.   verify_element_value(v,pattern*pattern);
  172.   v.sqrt();
  173.   verify_element_value(v,pattern);
  174.   v1.abs();
  175.   verify_element_value(v1,pattern);
  176.   assert( v == v1 );
  177.  
  178.   {
  179.     cout << "\nCheck out to see that sin^2(x) + cos^2(x) = 1\n";
  180.     for(i=v.q_lwb(); i<=v.q_upb(); i++)
  181.       v(i) = 2*M_PI/v.q_no_elems() * i;
  182.     class SinAction : public ElementPrimAction
  183.     {
  184.       void operation(REAL& element) { element = sin(element); }
  185.       public: SinAction(void) {}
  186.     };
  187.     v.apply(SinAction());
  188.     Vector v2 = v;
  189.     
  190.     class CosAction : public ElementAction
  191.     {
  192.       double factor;
  193.       void operation(REAL& element) { element = cos(factor*i); }
  194.       public:
  195.       CosAction(const int no_elems): factor(2*M_PI/no_elems) {}
  196.     };
  197.     v1.apply(CosAction(v.q_no_elems()));
  198.     Vector v3 = v1;
  199.     v.sqr();
  200.     v1.sqr();
  201.     v += v1;
  202.     verify_element_value(v,1);
  203.  
  204.     cout << "\n\tdo it again through LazyMatrix promise of a vector" << endl;
  205.     class square_add : public LazyMatrix, public ElementAction
  206.     {
  207.       const Vector &v1; Vector &v2;
  208.       void operation(REAL& element)
  209.               { assert(j==1); element = v1(i)*v1(i) + v2(i)*v2(i); }
  210.      void fill_in(Matrix& m) const { m.apply(*this); }
  211.      public: square_add(const Vector& _v1, Vector& _v2) :
  212.        LazyMatrix(_v1.q_row_lwb(),_v1.q_row_upb(),1,1),
  213.        v1(_v1), v2(_v2) {}
  214.     };
  215.     Vector vres = square_add(v2,v3);
  216.     Vector vres1 = v2; assert( !(vres1 == vres) );
  217.     verify_element_value(vres,1);
  218.     vres1 = square_add(v2,v3);
  219.     verify_element_value(vres1,1);
  220.   }
  221.  
  222.   cout << "\nVerify constructor with initialization\n";
  223.   Vector vi(1,5,0.0,1.0,2.0,3.0,4.0,"END");
  224.   Vector vit(5);
  225.   for(i=vit.q_lwb(); i<=vit.q_upb(); i++)
  226.     vit(i) = i-1;
  227.   verify_vector_identity(vi,vit);
  228.   
  229.   cout << "\nDone\n\n";
  230. }
  231.  
  232. /*
  233.  *------------------------------------------------------------------------
  234.  *            Test binary vector operations
  235.  */
  236.  
  237. static void test_binary_op(const int vsize)
  238. {
  239.   const double pattern = M_PI;
  240.   register int i;
  241.  
  242.   cout << "\n---> Test Binary Vector operations\n";
  243.  
  244.   Vector v(2,vsize+1);
  245.   Vector v1(v);
  246.   Vector vp(v);
  247.  
  248.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  249.     vp(i) = (i-v.q_no_elems()/2.)*pattern;
  250.   
  251.  
  252.   cout << "\nVerify assignment of a vector to the vector\n";
  253.   v = pattern;
  254.   v1.clear();
  255.   v1 = v;
  256.   verify_element_value(v1,pattern);
  257.   assert( v1 == v );
  258.  
  259.   cout << "\nAdding one vector to itself, uniform pattern " << pattern << "\n";
  260.   v.clear(); v = pattern;
  261.   v1 = v; v1 += v1;
  262.   verify_element_value(v1,2*pattern);
  263.   cout << "  subtracting two vectors ...\n";
  264.   v1 -= v;
  265.   verify_element_value(v1,pattern);
  266.   cout << "  subtracting the vector from itself\n";
  267.   v1 -= v1;
  268.   verify_element_value(v1,0);
  269.   cout << "  adding two vectors together\n";
  270.   v1 += v;
  271.   verify_element_value(v1,pattern);
  272.  
  273.   cout << "\nArithmetic operations on vectors with not the same elements\n";
  274.   cout << "   adding vp to the zero vector...\n";
  275.   v.clear(); v += vp;
  276. //  verify_vector_identity(v,vp);
  277.   assert( v == vp );
  278.   v1 = v;
  279.   cout << "   making v = 3*vp and v1 = 3*vp, via add() and succesive mult\n";
  280.   add(v,2,vp);
  281.   v1 += v1; v1 += vp;
  282.   verify_vector_identity(v,v1);
  283.   cout << "   clear both v and v1, by subtracting from itself and via add()\n";
  284.   v1 -= v1;
  285.   add(v,-3,vp);
  286.   verify_vector_identity(v,v1);
  287.  
  288.   cout << "\nTesting element-by-element multiplications and divisions\n";
  289.   cout << "   squaring each element with sqr() and via multiplication\n";
  290.   v = vp; v1 = vp;
  291.   v.sqr();
  292.   element_mult(v1,v1);
  293.   verify_vector_identity(v,v1);
  294.   cout << "   compare (v = pattern^2)/pattern with pattern\n";
  295.   v = pattern; v1 = pattern;
  296.   v.sqr();
  297.   element_div(v,v1);
  298.   verify_element_value(v,pattern);
  299.   compare(v1,v,"Original vector and vector after squaring and dividing");
  300.  
  301.   cout << "\nDone\n";
  302. }
  303.  
  304. /*
  305.  *------------------------------------------------------------------------
  306.  *            Verify the norm calculation
  307.  */
  308.  
  309. static void test_norms(const int vsize)
  310. {
  311.   cout << "\n---> Verify norm calculations\n";
  312.  
  313.   const double pattern = 10.25;
  314.  
  315.   if( vsize % 2 == 1 )
  316.     _error("Sorry, size of the vector to test must be even for this test\n");
  317.  
  318.   Vector v(vsize);
  319.   Vector v1(v);
  320.  
  321.   cout << "\nAssign " << pattern << " to all the elements and check norms\n";
  322.   v = pattern;
  323.   cout << "  1. norm should be pattern*no_elems\n";
  324.   assert( v.norm_1() == pattern*v.q_no_elems() );
  325.   cout << "  Square of the 2. norm has got to be pattern^2 * no_elems\n";
  326.   assert( v.norm_2_sqr() == sqr(pattern)*v.q_no_elems() );
  327.   cout << "  Inf norm should be pattern itself\n";
  328.   assert( v.norm_inf() == pattern );
  329.   cout << "  Scalar product of vector by itself is the sqr(2. vector norm)\n";
  330.   assert( v.norm_2_sqr() == v * v );
  331.  
  332.   double ap_step = 1;
  333.   double ap_a0   = -pattern;
  334.   int n = v.q_no_elems();
  335.   cout << "\nAssign the arithm progression with 1. term " << ap_a0 <<
  336.           "\nand the difference " << ap_step << "\n";
  337.   register int i;
  338.   for(i=v.q_lwb(); i<=v.q_upb(); i++)
  339.     v(i) = (i-v.q_lwb())*ap_step + ap_a0;
  340.   int l = min(max((int)ceil(-ap_a0/ap_step),0),n);
  341.   double norm = (2*ap_a0 + (l+n-1)*ap_step)/2*(n-l) +
  342.     (-2*ap_a0-(l-1)*ap_step)/2*l;
  343.   cout << "  1. norm should be " << norm << "\n";
  344.   assert( v.norm_1() == norm );
  345.   norm = n*( sqr(ap_a0) + ap_a0*ap_step*(n-1) + sqr(ap_step)*(n-1)*(2*n-1)/6);
  346.   cout << "  Square of the 2. norm has got to be "
  347.           "n*[ a0^2 + a0*q*(n-1) + q^2/6*(n-1)*(2n-1) ], or " << norm << "\n";
  348.   assert( v.norm_2_sqr() == norm );
  349.   norm = max(abs(v(v.q_lwb())),abs(v(v.q_upb())));
  350.   cout << "  Inf norm should be max(abs(a0),abs(a0+(n-1)*q)), ie " << norm <<
  351.           "\n";
  352.   assert( v.norm_inf() == norm );
  353.   cout << "  Scalar product of vector by itself is the sqr(2. vector norm)\n";
  354.   assert( v.norm_2_sqr() == v * v );
  355.  
  356.   v1.clear();
  357.   compare(v,v1,"Compare the vector v with a zero vector");
  358.  
  359.   cout << "\nConstruct v1 to be orthogonal to v as v(n), -v(n-1), v(n-2)...\n";
  360.   for(i=0; i<v1.q_no_elems(); i++)
  361.     v1(i+v1.q_lwb()) = v(v.q_upb()-i) * ( i % 2 == 1 ? -1 : 1 );
  362.   cout << "||v1|| has got to be equal ||v|| regardless of the norm def\n";
  363.   assert( v1.norm_1() == v.norm_1() );
  364.   assert( v1.norm_2_sqr() == v.norm_2_sqr() );
  365.   assert( v1.norm_inf() == v.norm_inf() );
  366.   cout << "But the scalar product has to be zero\n";
  367.   assert( v1 * v == 0 );
  368.  
  369.   cout << "\nDone\n";
  370. }
  371.  
  372. /*
  373.  *------------------------------------------------------------------------
  374.  *        Test operations with vectors and matrix slices
  375.  */
  376.  
  377. static void test_matrix_slices(const int vsize)
  378. {
  379.   const REAL pattern = 8.625;
  380.   register int i;
  381.  
  382.   cout << "\n\n---> Test operations with vectors and matrix slices\n";
  383.  
  384.   Vector vc(0,vsize);
  385.   Vector vr(0,vsize+1);
  386.   Matrix m(0,vsize,0,vsize+1);
  387.  
  388.   cout << "\nCheck modifying the matrix column-by-column\n";
  389.   m = pattern;
  390.   assert( m == pattern );
  391.   for(i=m.q_col_lwb(); i <= m.q_col_upb(); i++)
  392.   {
  393.     MatrixColumn(m,i) = pattern-1;
  394.     assert( !( m == pattern ) && !( m != pattern ) );
  395.     MatrixColumn(m,i) *= 2;
  396.     vc = MatrixColumn(m,i);
  397.     verify_element_value(vc,2*(pattern-1));
  398.     vc = MatrixColumn(m, i+1 > m.q_col_upb() ? m.q_col_lwb() : i+1);
  399.     verify_element_value(vc,pattern);
  400.     MatrixColumn(m,i) *= 0.5;
  401.     MatrixColumn(m,i) += 1;
  402.     assert( m == pattern );
  403.   }
  404.  
  405.   assert( m == pattern );
  406.   for(i=m.q_col_lwb(); i <= m.q_col_upb(); i++)
  407.   {
  408.     vc = pattern+1;
  409.     MatrixColumn(m,i) = vc;
  410.     assert( !( m == pattern ) && !( m != pattern ) );
  411.     {
  412.       MatrixColumn mc(m,i);
  413.       for(register int j=m.q_row_lwb(); j<=m.q_row_upb(); j++)
  414.         mc(j) *= 4;
  415.     }
  416.     vc = MatrixColumn(m,i);
  417.     verify_element_value(vc,4*(pattern+1));
  418.     MatrixColumn(m,i) *= 0.25;
  419.     MatrixColumn(m,i) += -1;
  420.     vc = MatrixColumn(m,i);
  421.     verify_element_value(vc,pattern);
  422.     assert( m == pattern );
  423.   }
  424.  
  425.   cout << "\nCheck modifying the matrix row-by-row\n";
  426.   m = pattern;
  427.   assert( m == pattern );
  428.   for(i=m.q_row_lwb(); i <= m.q_row_upb(); i++)
  429.   {
  430.     MatrixRow(m,i) = pattern+2;
  431.     assert( !( m == pattern ) && !( m != pattern ) );
  432.     vr = MatrixRow(m,i);
  433.     verify_element_value(vr,pattern+2);
  434.     vr = MatrixRow(m,i+1 > m.q_row_upb() ? m.q_row_lwb() : i+1);
  435.     verify_element_value(vr,pattern);
  436.     MatrixRow(m,i) += -2;
  437.     assert( m == pattern );
  438.   }
  439.  
  440.   assert( m == pattern );
  441.   for(i=m.q_row_lwb(); i <= m.q_row_upb(); i++)
  442.   {
  443.     vr = pattern-2;
  444.     MatrixRow(m,i) = vr;
  445.     assert( !( m == pattern ) && !( m != pattern ) );
  446.     {
  447.       MatrixRow mr(m,i);
  448.       for(register int j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  449.         mr(j) *= 8;
  450.     }
  451.     vr = MatrixRow(m,i);
  452.     verify_element_value(vr,8*(pattern-2));
  453.     MatrixRow(m,i) *= 1./8;
  454.     MatrixRow(m,i) += 2;
  455.     vr = MatrixRow(m,i);
  456.     verify_element_value(vr,pattern);
  457.     assert( m == pattern );
  458.   }
  459.  
  460.   cout << "\nCheck modifying the matrix diagonal\n";
  461.   m = pattern;
  462.   (MatrixDiag)m = pattern-3;
  463.   assert( !( m == pattern ) && !( m != pattern ) );
  464.   vc = MatrixDiag(m);
  465.   verify_element_value(vc,pattern-3);
  466.   MatrixDiag(m) += 3;
  467.   assert( m == pattern );
  468.   vc = pattern+3;
  469.   (MatrixDiag)m = vc;
  470.   assert( !( m == pattern ) && !( m != pattern ) );
  471.   {
  472.     MatrixDiag md(m);
  473.     for(register int j=1; j<=md.q_ndiags(); j++)
  474.       md(j) /= 1.5;
  475.   }
  476.   vc = MatrixDiag(m);
  477.   verify_element_value(vc,(pattern+3)/1.5);
  478.   MatrixDiag(m) *= 1.5;
  479.   MatrixDiag(m) += -3;
  480.   vc = MatrixDiag(m);
  481.   verify_element_value(vc,pattern);
  482.   assert( m == pattern );
  483.  
  484.   cout << "\nCheck out to see that multiplying by diagonal is column-wise"
  485.           "\nmatrix multiplication\n";
  486.   Matrix mm(m);
  487.   Matrix m1(m.q_row_lwb(),::max(m.q_row_upb(),m.q_col_upb()),
  488.         m.q_col_lwb(),::max(m.q_row_upb(),m.q_col_upb()));
  489.   Vector vc1(vc), vc2(vc);
  490.   for(i=m.q_row_lwb(); i<m.q_row_upb(); i++)
  491.     MatrixRow(m,i) = pattern+i;        // Make a multiplicand
  492.   mm = m;                // Save it
  493.  
  494.   m1 = pattern+10;
  495.   for(i=vr.q_lwb(); i<=vr.q_upb(); i++)
  496.     vr(i) = i+2;
  497.   (MatrixDiag)m1 = vr;            // Make the other multiplicand
  498.   assert( !(m1 == pattern+10) );
  499.  
  500.   m *= MatrixDiag(m1);
  501.   for(i=m.q_col_lwb(); i<=m.q_col_upb(); i++)
  502.   {
  503.     vc1 = MatrixColumn(mm,i);
  504.     vc1 *= vr(i);            // Do a column-wise multiplication
  505.     vc2 = MatrixColumn(m,i);
  506.     verify_vector_identity(vc1, vc2);
  507.   }
  508.  
  509.   cout << "\nDone\n";
  510. }
  511.  
  512. /*
  513.  *------------------------------------------------------------------------
  514.  *                Root module
  515.  */
  516.  
  517. main()
  518. {
  519.   cout << "\n\n" << _Minuses << 
  520.           "\n\t\tVerify Operations on Vectors\n";
  521.  
  522.   test_allocation();
  523.   test_element_op(20);
  524.   test_binary_op(20);
  525.   test_norms(20);
  526.   test_matrix_slices(20);
  527. }
  528.  
  529.